#include "Matrix.h"

void Matrix::allocate_memory() {
    if (rows <= 0 || cols <= 0) {
        data = nullptr;
        return;
    }
    data = new double* [rows];
    for (int i = 0; i < rows; ++i) {
        data[i] = new double[cols]();  // Initialize with zeros
    }
}

void Matrix::deallocate_memory() {
    if (data) {
        for (int i = 0; i < rows; ++i) {
            delete[] data[i];
        }
        delete[] data;
    }
    rows = 0;
    cols = 0;
    data = nullptr;
}

void Matrix::copy_from(const Matrix& other) {
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            data[i][j] = other.data[i][j];
        }
    }
}

Matrix::Matrix(int rows, int cols, double initialValue)
    : rows(rows), cols(cols) {
    allocate_memory();
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            data[i][j] = initialValue;
        }
    }
}

Matrix::Matrix(const Matrix& other)
    : rows(other.rows), cols(other.cols) {
    allocate_memory();
    copy_from(other);
}

Matrix::~Matrix() {
    deallocate_memory();
}

Matrix Matrix::operator+(const Matrix& other) const {
    if (rows != other.rows || cols != other.cols) {
        return Matrix();
    }
    Matrix result(rows, cols);
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            result.data[i][j] = data[i][j] + other.data[i][j];
        }
    }
    return result;
}

Matrix Matrix::operator-(const Matrix& other) const {
    if (rows != other.rows || cols != other.cols) {
        return Matrix();
    }
    Matrix result(rows, cols);
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            result.data[i][j] = data[i][j] - other.data[i][j];
        }
    }
    return result;
}

Matrix Matrix::operator*(const Matrix& other) const {
    if (cols != other.rows) {
        return Matrix();
    }
    Matrix result(rows, other.cols);
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < other.cols; ++j) {
            double sum = 0.0;
            for (int k = 0; k < cols; ++k) {
                sum += data[i][k] * other.data[k][j];
            }
            result.data[i][j] = sum;
        }
    }
    return result;
}

Matrix& Matrix::operator=(const Matrix& other) {
    if (this == &other) return *this;
    deallocate_memory();
    rows = other.rows;
    cols = other.cols;
    allocate_memory();
    copy_from(other);
    return *this;
}

Matrix Matrix::transpose() const {
    Matrix result(cols, rows);
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            result.data[j][i] = data[i][j];
        }
    }
    return result;
}

Matrix Matrix::inverse(const Matrix& mat) {
    int n = mat.rows;
    if (n != mat.cols) return Matrix();

    Matrix augmented(n, 2 * n);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            augmented(i, j) = mat(i, j);
        }
        augmented(i, i + n) = 1.0;
    }

    for (int i = 0; i < n; ++i) {
        int pivot = i;
        for (int j = i; j < n; ++j) {
            if (std::abs(augmented(j, i)) > std::abs(augmented(pivot, i))) {
                pivot = j;
            }
        }

        if (std::abs(augmented(pivot, i)) < 1e-10) return Matrix();

        if (pivot != i) {
            for (int j = 0; j < 2 * n; ++j) {
                std::swap(augmented(i, j), augmented(pivot, j));
            }
        }

        double divisor = augmented(i, i);
        for (int j = 0; j < 2 * n; ++j) {
            augmented(i, j) /= divisor;
        }

        for (int k = 0; k < n; ++k) {
            if (k != i && std::abs(augmented(k, i)) > 1e-10) {
                double factor = augmented(k, i);
                for (int j = 0; j < 2 * n; ++j) {
                    augmented(k, j) -= factor * augmented(i, j);
                }
            }
        }
    }

    Matrix inverse(n, n);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            inverse(i, j) = augmented(i, j + n);
        }
    }
    return inverse;
}

double Matrix::operator()(int row, int col) const {
    if (row >= 0 && row < rows && col >= 0 && col < cols) {
        return data[row][col];
    }
    return 0.0;
}

double& Matrix::operator()(int row, int col) {
    static double dummy = 0.0;
    if (row >= 0 && row < rows && col >= 0 && col < cols) {
        return data[row][col];
    }
    return dummy;
}

void Matrix::print(const std::string& title) const {
    if (!title.empty()) {
        std::cout << "\n" << title << " (" << rows << "x" << cols << "):\n";
    }
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            std::cout << std::setw(8) << std::fixed << std::setprecision(2)
                << data[i][j];
        }
        std::cout << std::endl;
    }
}